/*********************************************************************/
/* Multiple Precision float (MPF) Fortran -> C API                   */
/*********************************************************************/
#include <stdlib.h>
#include <string.h>
#include "mpf.h"
#include "cfort.h"

#define OPT_CHUNK 10

static size_t opt_sz=0;
static int opt_indx=-1;
static mp_size_t __mpf_prec=MPF_BITS_TO_PREC (53);
static mpf_ptr **x=0;
static int *x_sz=0;
static mpf_ptr buf_mpf;
static mp_ptr tp;

static void __mpf_alloc (mpf_ptr *);
static void __mpf_free (mpf_ptr *);

#define mpf_set_prec FortranName(MPF_SET_PREC,mpf_set_prec)
extern void mpf_set_prec (int *, int *);

#define mpf_alloc FortranName(MPF_ALLOC,mpf_alloc)
extern void mpf_alloc (int *, int *, int *);

#define mpf_realloc FortranName(MPF_REALLOC,mpf_realloc)
extern void mpf_realloc (int *, int *, int *);

#define mpf_free FortranName(MPF_FREE,mpf_free)
extern void mpf_free (int *);

#define mpf_size FortranName(MPF_SIZE,mpf_size)
extern void mpf_size (int *, int *);

#define mpf_set FortranName(MPF_SET,mpf_set)
extern void mpf_set (int *, Real *, int *, int *, int *, int *);

#define mpf_setv FortranName(MPF_SETV,mpf_setv)
extern void mpf_setv (int *, Real *, int *, int *, int *, int *);

#define mpf_get FortranName(MPF_GET,mpf_get)
extern void mpf_get (int *, Real *, int *, int *, int *, int *);

#define mpf_getv FortranName(MPF_GETV,mpf_getv)
extern void mpf_getv (int *, Real *, int *, int *, int *, int *);

#define mpf_setraw FortranName(MPF_SETRAW,mpf_setraw)
extern void mpf_setraw (int *, int *, int *, unsigned char *, int *, int *);

#define mpf_getraw FortranName(MPF_GETRAW,mpf_getraw)
extern void mpf_getraw (int *, int *, int *, unsigned char *, int *, int *);

#define mpf_add FortranName(MPF_ADD,mpf_add)
extern void mpf_add (int *, Real *, int *, int *, int *, int *);

#define mpf_sub FortranName(MPF_SUB,mpf_sub)
extern void mpf_sub (int *, Real *, int *, int *, int *, int *);

#define mpf_addraw FortranName(MPF_ADDRAW,mpf_addraw)
extern void mpf_addraw (int *, int *, int *, unsigned char *, int *, int *);

void mpf_set_prec (int *prec_in_bit, int *impfs)
{
   __mpf_prec = MPF_BITS_TO_PREC ((unsigned long)*prec_in_bit);
   __mpf_alloc(&buf_mpf);
   tp = (mp_ptr)malloc ((__mpf_prec + 1) * BYTES_PER_MP_LIMB);
   *impfs = ((__mpf_prec + 1) * BYTES_PER_MP_LIMB) + 2*sizeof(int);
}

void mpf_alloc (int* dim_f, int *ihandle, int *ier)
{
   int i,l;

   *ier = -1;
   l = *dim_f;

   if (l==0) {
     *ihandle = 0;
     *ier = 0;
     return;
   }

   if (opt_indx+1==opt_sz)
     {
       opt_sz += OPT_CHUNK;
       x = (mpf_ptr **)realloc (x,opt_sz*sizeof(mpf_ptr *));
       if (!x) return;
       x_sz = (int *)realloc (x_sz,opt_sz*sizeof(int));
       if (!x_sz) return;
     }
   opt_indx ++;
   x[opt_indx] = (mpf_ptr *)malloc (l*sizeof(mpf_ptr));
   if (!x[opt_indx]) return;
   x_sz[opt_indx] = l;

   for (i=0;i<l;i++) {
     __mpf_alloc(&x[opt_indx][i]);
     if (!x[opt_indx][i]) return;
   }

   *ier = 0;
   *ihandle = opt_indx+1;
}

void mpf_realloc (int* dim_f, int *ihandle, int *ier)
{
   int i,l,indx;

   if (*ihandle==0) {
     mpf_alloc(dim_f,ihandle,ier);
     return;
   }

   *ier = -1;
   indx = *ihandle-1;
   l = *dim_f;

   if (l>x_sz[indx]) {
     x[indx] = (mpf_ptr *)realloc (x[indx],l*sizeof(mpf_ptr));
     if (!x[indx]) return;
     for (i=x_sz[indx];i<l;i++) {
       __mpf_alloc(&x[indx][i]);
       if (!x[indx][i]) return;
     }
   }
   else {
     for (i=l;i<x_sz[indx];i++)
       __mpf_free(&x[indx][i]);
   }
   
   x_sz[indx] = l;
   *ier = 0;
}

void mpf_free (int *ihandle)
{
   int i,l,indx;

   if (*ihandle==0)
     return;

   indx = *ihandle-1;
   l = x_sz[indx];

   for (i=0;i<l;i++)
     __mpf_free(&x[indx][i]);

   x_sz[indx] = 0;

}

void mpf_size (int *ihandle, int *size)
{
   int indx;

   *size = 0;
   if (*ihandle==0)
     return;

   indx = *ihandle-1;
   *size = x_sz[indx];

}

void mpf_set (int *n, Real *y, int *incy, int *xaddr, int *incx, int *ihandle)
{
   int i,j,k,m,indx;
   int negative;
   double d;
   mpf_ptr xp;

   j = *incy;
   indx = *ihandle-1;
   k = *xaddr-1;
   m = *incx;

   for (i=0;i<*n;i++)
     {
       d = y[i*j];
       negative = d < 0;
       xp = x[indx][k];
       __D2MPF(xp,ABS(d),negative);
       k += m;
     }
}

void mpf_setv (int *n, Real *y, int *incy, int *xaddr, int *incx, int *ihandle)
{
   int i,j,k,indx;
   int negative;
   double d;
   mpf_ptr xp;

   indx = *ihandle-1;
   k = *xaddr-1;
   j = 0;

   for (i=0;i<*n;i++)
     {
       d = y[j];
       negative = d < 0;
       xp = x[indx][k];
       __D2MPF(xp,ABS(d),negative);
       k += incx[i];
       j += incy[i];
     }
}

void mpf_get (int *n, Real *y, int *incy, int *xaddr, int *incx, int *ihandle)
{
   int i,j,k,m,indx;
   mp_size_t size,ii,n_limbs_to_use;
   int negative;
   mp_ptr qp;
   mp_exp_t exp;
   union ieee_double_extract ieee_x;
   mpf_ptr xp;
   double d;

   j = *incy;
   indx = *ihandle-1;
   k = *xaddr-1;
   m = *incx;

   for (i=0;i<*n;i++) {
       
     xp = x[indx][k];
     size = SIZ(xp);
     d = 0;

     if (size != 0) {

       exp = EXP(xp);
       qp = PTR(xp);
       negative = size < 0;
       size = ABS (size);
       ieee_x.d = qp[size - 1];
       n_limbs_to_use = MIN (LIMBS_PER_DOUBLE, size);
       for (ii = 2; ii <= n_limbs_to_use; ii++) {
         ieee_x.d = ieee_x.d * MP_BASE_AS_DOUBLE + qp[size - ii];
       }

       exp = (exp - n_limbs_to_use) * BITS_PER_MP_LIMB ;

#if ( _MPF_IEEE_FLOATS == 1 )

       exp += (mp_exp_t)ieee_x.s.exp ;
       ieee_x.s.exp = exp;

       if (exp >= 2047) {
         /* Return +-infinity */
         ieee_x.s.exp = 2047;
         ieee_x.s.manl = 0;
         ieee_x.s.manh = 0;
       }
       else if (exp < 1) {
         /* here we choose to flush denormalized values 
         to zero */
         ieee_x.d = 0;
       }
#else
       {

         double factor,r;

         factor = 2.0;

         if (exp < 0) {
           factor = 0.5;
           exp = -exp;
         }
         r = ieee_x.d;
         if (exp != 0) {
           if ((exp & 1) != 0) {
             r *= factor;
           }
           exp >>= 1;
           while (exp != 0) {
             factor *= factor;
             if ((exp & 1) != 0) {
               r *= factor;
             }
             exp >>= 1;
           }
         }
         ieee_x.d = r;
       }  
#endif
       d = negative ? -ieee_x.d : ieee_x.d;
     }

     y[i*j] = d;
     k += m;

   }
}

void mpf_getv (int *n, Real *y, int *incy, int *xaddr, int *incx, int *ihandle)
{
   int i,j,k,indx;
   mp_size_t size,ii,n_limbs_to_use;
   int negative;
   mp_ptr qp;
   mp_exp_t exp;
   union ieee_double_extract ieee_x;
   mpf_ptr xp;
   double d;

   indx = *ihandle-1;
   k = *xaddr-1;
   j = 0;

   for (i=0;i<*n;i++) {
       
     xp = x[indx][k];
     size = SIZ(xp);
     d = 0;

     if (size != 0) {

       exp = EXP(xp);
       qp = PTR(xp);
       negative = size < 0;
       size = ABS (size);
       ieee_x.d = qp[size - 1];
       n_limbs_to_use = MIN (LIMBS_PER_DOUBLE, size);
       for (ii = 2; ii <= n_limbs_to_use; ii++) {
         ieee_x.d = ieee_x.d * MP_BASE_AS_DOUBLE + qp[size - ii];
       }

       exp = (exp - n_limbs_to_use) * BITS_PER_MP_LIMB ;

#if ( _MPF_IEEE_FLOATS == 1 )

       exp += (mp_exp_t)ieee_x.s.exp ;
       ieee_x.s.exp = exp;

       if (exp >= 2047) {
         /* Return +-infinity */
         ieee_x.s.exp = 2047;
         ieee_x.s.manl = 0;
         ieee_x.s.manh = 0;
       }
       else if (exp < 1) {
         /* here we choose to flush denormalized values 
         to zero */
         ieee_x.d = 0;
       }
#else
       {

         double factor,r;

         factor = 2.0;

         if (exp < 0) {
           factor = 0.5;
           exp = -exp;
         }
         r = ieee_x.d;
         if (exp != 0) {
           if ((exp & 1) != 0) {
             r *= factor;
           }
           exp >>= 1;
           while (exp != 0) {
             factor *= factor;
             if ((exp & 1) != 0) {
               r *= factor;
             }
             exp >>= 1;
           }
         }
         ieee_x.d = r;
       }  
#endif
       d = negative ? -ieee_x.d : ieee_x.d;
     }

     y[j] = d;
     k += incx[i];
     j += incy[i];

   }
}

void mpf_add (int *n, Real *y, int *incy, int *xaddr, int *incx, int *ihandle)
{
   int       i,j,k,m,indx;
   int       negative;
   double    bufd;
   mp_size_t rsize,usize,size;
   mp_size_t ediff;
   mp_exp_t  rexp,uexp;
   mp_ptr    rp,up;
   mpf_ptr   xp;

   j = *incy;
   indx = *ihandle-1;
   k = *xaddr-1;
   m = *incx;

   for (i=0;i<*n;i++)
     {
       bufd = y[i*j];
       negative = bufd < 0;
       __D2MPF(buf_mpf,ABS(bufd),negative);

       xp = x[indx][k];

       usize = SIZ(buf_mpf);
       if (usize==0)
         goto done;

       rsize = SIZ(xp);
       if (rsize==0)
         {
           xp = buf_mpf; buf_mpf = x[indx][k];
           goto done;
         }

       rexp = EXP(xp);
       uexp = EXP(buf_mpf);

       /* If signs of operands are different, perform substraction.  */
       if ((rsize^usize) < 0)
         {
           /* Signs are now known to be the same.  */
           negative = rsize < 0;

           /* Make r to be the operand with the largest exponent.  */
           if (rexp < uexp)
             {
                mpf_ptr t;
                mp_exp_t te;
                mp_size_t ts;
                t = xp; xp = buf_mpf; buf_mpf = t;
                te = rexp; rexp = uexp; uexp = te;
                ts = rsize; rsize = usize; usize = ts;
                negative ^= 1;
             }

           rsize = ABS(rsize);
           usize = ABS(usize);

           up = PTR(buf_mpf);
           rp = PTR(xp);
           ediff = rexp - uexp;

           __SUB;

         }
       else /* performing addition */
         {
           /* Signs are now known to be the same.  */
           negative = rsize < 0;

           /* Make r to be the operand with the largest exponent.  */
           if (rexp < uexp)
             {
                mpf_ptr t;
                mp_exp_t te;
                mp_size_t ts;
                t = xp; xp = buf_mpf; buf_mpf = t;
                te = rexp; rexp = uexp; uexp = te;
                ts = rsize; rsize = usize; usize = ts;
             }

           rsize = ABS(rsize);
           usize = ABS(usize);

           up = PTR(buf_mpf);
           rp = PTR(xp);
           ediff = rexp - uexp;

           __ADD;

         }
       done:
       x[indx][k] = xp;
       k += m;
     }
}

void mpf_sub (int *n, Real *y, int *incy, int *xaddr, int *incx, int *ihandle)
{
   int       i,j,k,m,indx;
   int       negative;
   double    bufd;
   mp_size_t rsize,usize,size;
   mp_size_t ediff;
   mp_exp_t  rexp,uexp;
   mp_ptr    rp,up;
   mpf_ptr   xp;

   j = *incy;
   indx = *ihandle-1;
   k = *xaddr-1;
   m = *incx;

   for (i=0;i<*n;i++)
     {
       bufd = y[i*j];
       negative = bufd < 0;
       __D2MPF(buf_mpf,ABS(bufd),negative);

       xp = x[indx][k];

       usize = SIZ(buf_mpf);
       if (usize==0)
         goto done;

       rsize = SIZ(xp);
       if (rsize==0)
         {
           SIZ(buf_mpf) = -usize;
           xp = buf_mpf; buf_mpf = x[indx][k];
           goto done;
         }

       rexp = EXP(xp);
       uexp = EXP(buf_mpf);

       /* If signs of operands are the same, perform substraction. */
       if ((rsize^usize) >= 0)
         {
           /* Signs are now known to be the same.  */
           negative = rsize < 0;

           /* Make r to be the operand with the largest exponent.  */
           if (rexp < uexp)
             {
                mpf_ptr t;
                mp_exp_t te;
                mp_size_t ts;
                t = xp; xp = buf_mpf; buf_mpf = t;
                te = rexp; rexp = uexp; uexp = te;
                ts = rsize; rsize = usize; usize = ts;
                negative ^= 1;
             }

           rsize = ABS(rsize);
           usize = ABS(usize);

           up = PTR(buf_mpf);
           rp = PTR(xp);
           ediff = rexp - uexp;

           __SUB;

         }
       else /* performing addition */
         {
           /* Signs are now known to be the same.  */
           negative = rsize < 0;

           /* Make r to be the operand with the largest exponent.  */
           if (rexp < uexp)
             {
                mpf_ptr t;
                mp_exp_t te;
                mp_size_t ts;
                t = xp; xp = buf_mpf; buf_mpf = t;
                te = rexp; rexp = uexp; uexp = te;
                ts = rsize; rsize = usize; usize = ts;
             }

           rsize = ABS(rsize);
           usize = ABS(usize);

           up = PTR(buf_mpf);
           rp = PTR(xp);
           ediff = rexp - uexp;

           __ADD;

         }
       done:
       x[indx][k] = xp;
       k += m;
     }
}


void mpf_setraw (int *n, int *xaddr, int *incx, unsigned char *stream, int *lstream, int *ihandle)
{
   int i,j,k,m,indx;
   mp_size_t s;
   mp_size_t rsize;
   mp_ptr rp;
   mpf_ptr xp;
   int out_bytesize,out_exp;
   mp_limb_t hi_limb;
   int negative;
   int n_bytes_in_hi_limb;
   unsigned char *cp = stream + *lstream;

   indx = *ihandle-1;
   k = *xaddr-1;
   m = *incx;

   for (i=0;i<*n;i++)
     {
       xp = x[indx][k];
       rsize = SIZ(xp);
       if (rsize == 0)
         {
           for (j = 4 - 1; j >= 0; j--)
             *cp++ = 0;
         }
       else
         {
           negative = rsize < 0;
           rsize = ABS(rsize);
           out_exp = EXP(xp);
           rp = PTR(xp);
           hi_limb = rp[rsize - 1];
           for (j = BYTES_PER_MP_LIMB - 1; j > 0; j--)
             {
               if ((hi_limb >> j * BITS_PER_CHAR) != 0)
                 break;
             }
           n_bytes_in_hi_limb = j + 1;
           out_bytesize = BYTES_PER_MP_LIMB * (rsize - 1) + n_bytes_in_hi_limb;

           out_bytesize = negative ? -out_bytesize : out_bytesize;
           /* Make the size 4 bytes on all machines, to make the format portable. */
           for (j = 4 - 1; j >= 0; j--)
             *cp++ = (out_bytesize >> (j * BITS_PER_CHAR)) % (1 << BITS_PER_CHAR) ;
           /* Make the size 4 bytes on all machines, to make the format portable. */
           for (j = 4 - 1; j >= 0; j--)
             *cp++ = (out_exp >> (j * BITS_PER_CHAR)) % (1 << BITS_PER_CHAR) ;

           /* Output from the most significant limb to the least significant limb,
              with each limb also output in decreasing significance order.  */

           /* Output the most significant limb separately, since we will only
              output some of its bytes.  */
           for (j = n_bytes_in_hi_limb - 1; j >= 0; j--)
             *cp++ = (hi_limb >> (j * BITS_PER_CHAR)) % ((mp_limb_t)1 << BITS_PER_CHAR);

           /* Output the remaining limbs. */
           for (s = rsize - 2; s >= 0; s--)
             {
               mp_limb_t x_limb;

               x_limb = rp[s];
               for (j = BYTES_PER_MP_LIMB - 1; j >= 0; j--)
                 *cp++ = (x_limb >> (j * BITS_PER_CHAR)) % ((mp_limb_t)1 << BITS_PER_CHAR);
             }
         }
       k += m;
     }

     *lstream = (int)(cp - stream);
}

void mpf_getraw (int *n, int *xaddr, int *incx, unsigned char *stream, int *lstream, int *ihandle)
{
   int i,j,k,m,indx;
   mp_size_t s;
   mp_size_t rsize;
   mpf_ptr xp;
   mp_ptr rp;
   unsigned int c;
   mp_limb_t x_limb;
   int in_bytesize,in_exp;
   int negative;
   unsigned char *cp = stream + *lstream;

   indx = *ihandle-1;
   k = *xaddr-1;
   m = *incx;

   for (i=0;i<*n;i++)
     {
       xp = x[indx][k];
       in_bytesize = 0;
       /* Read 4-byte size */
       for (j = 4 - 1; j >= 0; j--)
         {
           c = *cp++;
           in_bytesize = (in_bytesize << BITS_PER_CHAR) | c;
         }
       negative = in_bytesize < 0;
       in_bytesize = ABS (in_bytesize);
       rsize = (in_bytesize + BYTES_PER_MP_LIMB - 1) / BYTES_PER_MP_LIMB;
       if (rsize == 0)
         {
           SIZ(xp) = 0;
         }
       else
         {
           rp = PTR(xp);
           /* Read 4-byte exponent */
           in_exp = 0;
           for (j = 4 - 1; j >= 0; j--)
             {
               c = *cp++;
               in_exp = (in_exp << BITS_PER_CHAR) | c;
             }
           EXP(xp) = in_exp;

           x_limb = 0;
           for (j = (in_bytesize - 1) % BYTES_PER_MP_LIMB; j >= 0; j--)
             {
               c = *cp++;
               x_limb = (x_limb << BITS_PER_CHAR) | c;
             }
           rp[rsize - 1] = x_limb;

           for (s = rsize - 2; s >= 0; s--)
             {
               x_limb = 0;
               for (j = BYTES_PER_MP_LIMB - 1; j >= 0; j--)
                 {
                   c = *cp++;
                   x_limb = (x_limb << BITS_PER_CHAR) | c;
                 }
               rp[s] = x_limb;
             }

           SIZ(xp) = negative ? -rsize : rsize;
         }
       k += m;
     }

     *lstream = (int)(cp - stream);
}

void mpf_addraw (int *n, int *xaddr, int *incx, unsigned char *stream, int *lstream, int *ihandle)
{
   int i,j,k,m,indx;
   int negative;
   mp_size_t rsize,usize,size;
   mp_size_t ediff;
   mp_exp_t  rexp,uexp;
   mp_ptr    rp,up;
   mp_limb_t x_limb;
   mpf_ptr xp;
   mp_size_t s;
   unsigned int c;
   int in_bytesize,in_exp;
   unsigned char *cp = stream + *lstream;

   indx = *ihandle-1;
   k = *xaddr-1;
   m = *incx;

   for (i=0;i<*n;i++)
     {
       xp = x[indx][k];

       in_bytesize = 0;
       /* Read 4-byte size */
       for (j = 4 - 1; j >= 0; j--)
         {
           c = *cp++;
           in_bytesize = (in_bytesize << BITS_PER_CHAR) | c;
         }
       negative = in_bytesize < 0;
       in_bytesize = ABS (in_bytesize);
       usize = (in_bytesize + BYTES_PER_MP_LIMB - 1) / BYTES_PER_MP_LIMB;
       if (usize == 0) goto done;

       up = PTR(buf_mpf);
       /* Read 4-byte exponent */
       in_exp = 0;
       for (j = 4 - 1; j >= 0; j--)
         {
           c = *cp++;
           in_exp = (in_exp << BITS_PER_CHAR) | c;
         }
       EXP(buf_mpf) = in_exp;

       x_limb = 0;
       for (j = (in_bytesize - 1) % BYTES_PER_MP_LIMB; j >= 0; j--)
         {
           c = *cp++;
           x_limb = (x_limb << BITS_PER_CHAR) | c;
         }
       up[usize - 1] = x_limb;

       for (s = usize - 2; s >= 0; s--)
         {
           x_limb = 0;
           for (j = BYTES_PER_MP_LIMB - 1; j >= 0; j--)
             {
               c = *cp++;
               x_limb = (x_limb << BITS_PER_CHAR) | c;
             }
           up[s] = x_limb;
         }

       usize = negative ? -usize : usize;
       SIZ(buf_mpf) = usize;

       rsize = SIZ(xp);
       if (rsize==0)
         {
           xp = buf_mpf; buf_mpf = x[indx][k];
           goto done;
         }

       /* If signs of operands are different, perform substraction.  */
       if ((rsize^usize) < 0)
         {
           /* Signs are now known to be the same.  */
           negative = rsize < 0;
           rexp = EXP(xp);
           uexp = EXP(buf_mpf);

           /* Make r to be the operand with the largest exponent.  */
           if (rexp < uexp)
             {
                mpf_ptr t;
                mp_exp_t te;
                mp_size_t ts;
                t = xp; xp = buf_mpf; buf_mpf = t;
                te = rexp; rexp = uexp; uexp = te;
                ts = rsize; rsize = usize; usize = ts;
                negative ^= 1;
             }

           rsize = ABS(rsize);
           usize = ABS(usize);

           up = PTR(buf_mpf);
           rp = PTR(xp);
           ediff = rexp - uexp;

           __SUB;

         }
       else /* performing addition */
         {
           /* Signs are now known to be the same.  */
           negative = rsize < 0;
           rexp = EXP(xp);
           uexp = EXP(buf_mpf);

           /* Make r to be the operand with the largest exponent.  */
           if (rexp < uexp)
             {
                mpf_ptr t;
                mp_exp_t te;
                mp_size_t ts;
                t = xp; xp = buf_mpf; buf_mpf = t;
                te = rexp; rexp = uexp; uexp = te;
                ts = rsize; rsize = usize; usize = ts;
             }

           rsize = ABS(rsize);
           usize = ABS(usize);

           up = PTR(buf_mpf);
           rp = PTR(xp);
           ediff = rexp - uexp;

           __ADD;

         }
       done:
       x[indx][k] = xp;
       k += m;
     }

     *lstream = (int)(cp - stream);
}

/*****************************************************************************/
/* Internal services                                                         */
/*****************************************************************************/

static void __mpf_alloc(mpf_ptr *r)
{
   *r = (mpf_ptr)malloc (sizeof(MP_FLOAT));
   if (!*r) return;
   PTR(*r)  = (mp_ptr)malloc ((__mpf_prec + 1) * BYTES_PER_MP_LIMB);
   PREC(*r) = __mpf_prec;
   SIZ(*r)  = 0;
   EXP(*r)  = 0;
}

static void __mpf_free(mpf_ptr *r)
{
   free(PTR(*r));
   free(*r);
   *r = 0;
}
